--- title: Hotmax spectra keywords: fastai sidebar: home_sidebar summary: "Cherry picking the best spectra" description: "Cherry picking the best spectra" nb_path: "notebooks/40_hotmax.ipynb" ---
A MA-XRF spectral image data cube contains typically a million spectra or so. Far too many to analyze them individually. Furthermore, most of these spectra have low signals and are rather noisy. An approach to overcome these computational problems is to base a further data analysis solely on the max spectrum. Such an approach is used e.g. by DataMuncher developed by Matthias Alfeld. As discussed above, the max spectrum essentially is an envelope function that provides a highly informative summary of all spectra. Different peaks in the max spectrum envelope can originate from from different spatial positions in the spectral image. Although this mathematical combination of spectral signals reduces computational analysis time, down the road it complicates our task of unmixing spectra and attribution of specific peak patterns to different chemical elements.
A more sophisticated approach is to locate for each peak in the max spectrum envelope which specific pixel spectrum is responsible for that specific peak. Loosely speaking, which pixels in the spectral image data cube are 'hotmax'? Another way to explain this is to find the specific 'hotmax' spectra that touches the corresponding max spectrum peaks.
A requisite step in the data analysis now is to find the hotmax pixels and spectra. Locating them takes a few minutes, and should be done once using the get_hotmax_spectra() function. The user is prompted to inspect and save the result in the datastack file.
from maxrf4u import compute_hotmax_spectra
compute_hotmax_spectra('RP-T-1898-A-3689.datastack')
In further analysis our stored hotmax pixel indexes and hotmax spectra can now be accessed from file using the DataStack.read(<datapath>)) methods. Let's inspect the shapes of these arrays first.
from maxrf4u import DataStack
ds = DataStack('RP-T-1898-A-3689.datastack')
x_keVs = ds.read('maxrf_energies')
hotmax_spectra = ds.read('hotmax_spectra')
hotmax_spots = ds.read('hotmax_spots')
hotmax_pixels = ds.read('hotmax_pixels')
hotmax_mask = ds.read('hotmax_mask')
n_pixels, n_dims = hotmax_pixels.shape
n_spectra, n_channels = hotmax_spectra.shape
print(f'The array of hotmax pixels (actually voxels) contains {n_pixels} pixels.')
print(f'Each pixel consists of {n_dims} (y,x,z) indexes.\n')
print(f'There are {n_spectra} unique hotmax spots and corresponding spectra.')
print(f'Each spectrum contains {n_channels} intensity values.')
To get an idea, let's plot both the spatial (x, y) locations of the hotmax spots on the image of the drawing, and their corresponding hotmax spectra.
Now let's pick one of the 22 hotmax spectra and take a closer look. Each hotmax spectrum has one or more peaks (square red marker) that correspond to specific peaks in the the max spectrum envelope.
The measured spectrum consists of a a slowly varying baseline, the so-called the continuum ridge, with peaks added on top. Some of these peaks are caused by x-ray fluorescence of specific chemical elements or Other peaks are just noise. The difficulty with these kind of spectra is that the noise level varies with the signal according to Poisson statistics. The variance of noise is linearly proportional to the signal level. For this reason, we observe more noise on top of the continuum ridge.
In order to distinguish significant peaks from noise it is possible to estimate a baseline with a noise envelope for the hotmax spectra with the compute_hotmax_noise()function. The algorithm used in estimating the baseline function is the rolling ball filter. In accordance with Poisson statistics, the noise level is estimated as a square root function of the baseline level.
from maxrf4u import compute_hotmax_noise
compute_hotmax_noise('RP-T-1898-A-3689.datastack')
We can now create a HotmaxAtlas object to plot all spectra with their noise envelopes and sub peaks with the HotmaxAtlas.plot_spectra()) method.
from maxrf4u import HotmaxAtlas
hma = HotmaxAtlas('RP-T-1898-A-3689.datastack')
hma.plot_spectra()
Or we can plot a single spectrum with HotmaxAtlas.plot_spectrum(<nr>)) to inspect a specific hotmax spectrum.
hma.plot_spectrum(12);
Each specific chemical element present is known to cause a specific peak pattern. Now we need to solve the inverse problem, which I call 'the peak pattern puzzle'. Given a spectrum, which combination of elements is likely to have caused the combination of peaks we find in the spectrum? This is the topic of the next section...